Dominant Species Plot Editing

Testing alternative plots for dominant species responses

Code

# load packages
library(here)
library(readxl)
library(ggstream)
library(ggforce)
library(tidyverse)
library(gmRi)
library(scales)
library(rcartocolor)
library(patchwork)
library(gt)
library(ggpubr)
library(ggh4x)


# Set theme
theme_set(
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold",size = 14),
    axis.text = element_text(size = 11),
    axis.title = element_text(size = 12),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 11),
    strip.background = element_rect(color = "black"),
    panel.background = element_rect(color = "black"),
    panel.grid.minor = element_blank())
)

# levels for season factors
season_lvls <- c("Winter", "Spring", "Summer", "Autumn")


#
set.seed(123)

Code

gmRi::use_gmri_style_rmd()

Code

# Load the Data:
dat_path <- here("maria_data/maria_results_2023/")


# Function to load multiple files at once to a list
load_folder <- function(folder, col_names = TRUE){
  folder_path <- str_c(dat_path, folder)
  fnames <- list.files(folder_path, full.names = TRUE) %>% 
    setNames(str_remove_all(list.files(folder_path), ".csv")) 
  flist <- map(fnames, ~read_csv(.x, guess_max = 1e5, col_types = cols(), col_names = col_names))
  return(flist)
}
  
  

# Function to scale subgroups as a vector
scale_this <- function(x){ (x - mean(x, na.rm=TRUE)) / sd(x, na.rm=TRUE) }

Figure Edits for Figure 3 of the Simulation Study

Reviewer Comment: > In figure 3 its hard to see which temperature optima dominate, is there another way to visualize this pattern? Also using a common color scale for sizes of protists and copepods seems unnecessary, and reduces the dynamic range within each panel.

Fig 5. Dominant Species Biomass

Absolute seasonal biomass concentration anomalies of the dominant groups of (a.) protists, (b.) passive and (a.) active copepod feeders for the autumn HW scenario during the heatwave and 2 years after the heatwave.

Code

# | label: new-dom-species

# Maria res-supplied these tables for dominant species plots, these should fix
# the problems we were encountering with the reshaped data previously

dom_new_ac <- read_csv(here("maria_data/maria_results_2023/dom_taxa_reshaped/Active_Copepods_Dominant.csv")) %>% 
  select(1:9)
dom_new_pc <- read_csv(here("maria_data/maria_results_2023/dom_taxa_reshaped/Passive_Copepods_Dominant.csv"))
dom_new_pr <- read_csv(here("maria_data/maria_results_2023/dom_taxa_reshaped/Protists_Dominant.csv"))


# Just need to put the control data ahead of the seasons for each to capture starting conditions
# Might be tricky to make sure the taxa_id's go into the right ones

# So annoying dolyr doesn't want to support this function...
named_group_split <- function(.tbl, ...) {
  grouped <- group_by(.tbl, ...)
  names <- rlang::inject(paste(!!!group_keys(grouped), sep = " / "))

  grouped %>% 
    group_split() %>% 
    rlang::set_names(names)}



# -----------------

# Reassemble them with control out front
dom_new_ac <- dom_new_ac %>% 
  filter(hw_season != "Control") %>% 
  named_group_split(hw_season, season) %>% 
  map_dfr(function(x_grouping){
    
    # Get the details of the group we're appending onto:
    taxa_ids <- unique(x_grouping$taxa_id)
    hw_season_lab <- unique(x_grouping$hw_season)
    season_lab <- unique(x_grouping$season)
    
    # Pull the appropriate control, 
    # filter to the taxa in this group
    # label the hw_season and season so we can facet
    control_bio <- dom_new_ac %>% 
      filter(
        hw_season == "Control",
        taxa_id %in% taxa_ids) %>% 
      mutate(
        #season = season_lab,
        hw_season = hw_season_lab)
    
    # Append it back
    df_out <- bind_rows(control_bio, x_grouping)
    return(df_out)})  %>% 
  mutate(
    temp_opt = str_c(temp_opt, "C"),
    hw_season = factor(
      str_c(hw_season, " Heatwave"),
      levels = str_c(season_lvls, " Heatwave")),
    season = factor(season, levels = season_lvls))





# Reassemble them with control out front
dom_new_pc <- dom_new_pc %>% 
  filter(hw_season != "Control") %>% 
  named_group_split(hw_season, season) %>% 
  map_dfr(function(x_grouping){
    
    # Get the details of the group we're appending onto:
    taxa_ids <- unique(x_grouping$taxa_id)
    hw_season_lab <- unique(x_grouping$hw_season)
    season_lab <- unique(x_grouping$season)
    
    # Pull the appropriate control, 
    # filter to the taxa in this group
    # label the hw_season and season so we can facet
    control_bio <- dom_new_pc %>% 
      filter(
        hw_season == "Control",
        taxa_id %in% taxa_ids) %>% 
      mutate(
        #season = season_lab,
        hw_season = hw_season_lab)
    
    # Append it back
    df_out <- bind_rows(control_bio, x_grouping)
    return(df_out)})  %>% 
  mutate(
    temp_opt = str_c(temp_opt, "C"),
    hw_season = factor(
      str_c(hw_season, " Heatwave"),
      levels = str_c(season_lvls, " Heatwave")),
    season = factor(season, levels = season_lvls))





# Reassemble them with control out front
dom_new_pr <- dom_new_pr %>% 
  filter(hw_season != "Control") %>% 
  named_group_split(hw_season, season) %>% 
  map_dfr(function(x_grouping){
    
    # Get the details of the group we're appending onto:
    taxa_ids <- unique(x_grouping$taxa_id)
    hw_season_lab <- unique(x_grouping$hw_season)
    season_lab <- unique(x_grouping$season)
    
    # Pull the appropriate control, 
    # filter to the taxa in this group
    # label the hw_season and season so we can facet
    control_bio <- dom_new_pr %>% 
      filter(
        hw_season == "Control",
        taxa_id %in% taxa_ids) %>% 
      mutate(
        #season = season_lab,
        hw_season = hw_season_lab)
    
    # Append it back
    df_out <- bind_rows(control_bio, x_grouping)
    return(df_out)})  %>% 
  mutate(
    temp_opt = str_c(temp_opt, "C"),
    hw_season = factor(
      str_c(hw_season, " Heatwave"),
      levels = str_c(season_lvls, " Heatwave")),
    season = factor(season, levels = season_lvls))

Reviewer Response: 2 Color Scales

In figure 3 1. Its hard to see which temperature optima dominate, 2. is there another way to visualize this pattern? 3. Also using a common color scale for sizes of protists and copepods seems unnecessary, and reduces the dynamic range within each panel.

Code

# Table for hw rectangles
hw_fill <- c("transparent", "#FDDBC770")
rect_df <- data.frame(
  hw_status = c("Year 0:\nPre-heatwave\nEquilibrium", "Year 1:\nSeasonal heatwave"),
  xmin = c(-0.5, 0.5),
  xmax = c(0.5, 1.5),
  ymin = rep(-Inf, 2),
  ymax = rep(Inf, 2))

# shape assignment by temp
temp_shape_vals <- c("16C" = 15, "20C" = 16, "24C" = 17, "28C" = 18)


# Line Chart:


# Dominant Copepod
f3_a <- dom_new_ac %>% 
  ggplot() +
    geom_rect(
      data = rect_df,
      aes(xmin = xmin, xmax = xmax,
          ymin = ymin, ymax = ymax,
          fill = hw_status),
      color = "transparent",
      alpha = 0.9) +
    geom_line(
      aes(year, absolute_bio, group = taxa_id, color = max_size), 
      show.legend = F, 
      linewidth = 0.7, 
      alpha = 0.4) +
    geom_point(
      aes(year, absolute_bio, shape = fct_rev(as.character(temp_opt)), color = max_size), 
    size = 2) +
    scale_shape_manual(
      values = temp_shape_vals) +
    scale_color_carto_c(
      palette = "ag_Sunset",
      direction = -1,
      trans = "log10",
      labels = label_log(base = 10),
      breaks = 10^seq(-3, 3, 1),
      limits = 10^c(-2, 3),
      oob = oob_squish) +
    scale_fill_manual(values = hw_fill) +
    scale_x_continuous(
      limits = c(0, 7),
      expand = expansion(add = c(0.5,0.5)),
      breaks = c(0:9)) +
    guides(
      color = guide_colorbar(
        order = 1,
        title.position = "top", 
        title.hjust = 0, 
        direction = "horizontal",
        barwidth = unit(4, "cm"), 
        ticks.colour = "black", 
        frame.colour = "black"),
      shape = guide_legend(
        order = 2,
        nrow = 1,
        title.position = "top", 
        title.hjust = 0, 
        override.aes = list(size = 3)),
      fill = guide_legend(
        order = 3,
        title.position = "top", 
        title.hjust = 0, nrow = 1,
        keyheight = unit(1, "cm"),
        override.aes = list(color = "black"))) +
    facet_grid(season~hw_season, scales = "free_y") +
    theme(
      legend.position = "right",
      panel.grid = element_blank(),
      panel.background = element_rect(color = "black", fill = "transparent")) +
    labs(
      y = expression(paste("Absolute biomass (mg C ", m^{-3}, ")")),
      x = NULL,
      color = "Copepod Maximum Body Size",
      shape = "Temperature Optima",
      fill = "Heatwave Timing:"
    )




# Passive Copepod
f3_b <- dom_new_pc %>% 
  ggplot() +
    geom_rect(
      data = rect_df,
      aes(xmin = xmin, xmax = xmax,
          ymin = ymin, ymax = ymax,
          fill = hw_status),
      color = "transparent",
      alpha = 0.9) +
    geom_line(
      aes(year, absolute_bio, group = taxa_id, color = max_size), 
      show.legend = F, 
      linewidth = 0.7, 
      alpha = 0.4) +
    geom_point(
      aes(year, absolute_bio, shape = fct_rev(as.character(temp_opt)), color = max_size), 
    size = 2) +
    scale_shape_manual(
      values = temp_shape_vals) +
    scale_color_carto_c(
      palette = "ag_Sunset",
      direction = -1,
      trans = "log10",
      labels = label_log(base = 10),
      breaks = 10^seq(-3, 3, 1),
      limits = 10^c(-2, 3),
      oob = oob_squish) +
    scale_fill_manual(values = hw_fill) +
    scale_x_continuous(
      limits = c(0, 7),
      expand = expansion(add = c(0.5,0.5)),
      breaks = c(0:9)) +
    guides(
      shape = "none",
      color = "none",
      fill = "none") +
    facet_grid(season~hw_season, scales = "free_y") +
    theme(
      legend.position = "right",
      panel.grid = element_blank(),
      panel.background = element_rect(color = "black", fill = "transparent")) +
    labs(
      y = expression(paste("Absolute biomass (mg C ", m^{-3}, ")")),
      x = NULL,
      color = "Copepod Maximum Body Size",
      shape = "Temperature Optima",
      fill = "Heatwave Timing:"
    )




# Protists
f3_c <- dom_new_pr %>% 
  ggplot() +
    geom_rect(
      data = rect_df,
      aes(xmin = xmin, xmax = xmax,
          ymin = ymin, ymax = ymax,
          fill = hw_status),
      color = "transparent",
      alpha = 0.9) +
    geom_line(
      aes(year, absolute_bio, group = taxa_id, color = max_size), 
      show.legend = F, 
      linewidth = 0.7, 
      alpha = 0.4) +
    geom_point(
      aes(year, absolute_bio, shape = fct_rev(as.character(temp_opt)), color = max_size), 
      size = 2) +
    scale_shape_manual(
      values = temp_shape_vals) +
    scale_color_carto_c(
      palette = "ag_GrnYl",
      direction = -1,
      trans = "log10",
      labels = label_log(base = 10),
      breaks = 10^seq(-6, 3, 1),
      limits = 10^c(-6, -2),
      oob = oob_squish) +
    scale_fill_manual(values = hw_fill) +
    scale_x_continuous(
      limits = c(0, 7),
      expand = expansion(add = c(0.5,0.5)),
      breaks = c(0:9)) +
    guides(
      color = guide_colorbar(
        order = 1,
        title.position = "top", 
        title.hjust = 0, 
        direction = "horizontal",
        barwidth = unit(4, "cm"), 
        ticks.colour = "black", 
        frame.colour = "black"),
      shape = "none",
      fill = "none") +
    facet_grid(season~hw_season, scales = "free_y") +
    theme(
      legend.position = "right",
      panel.grid = element_blank(),
      panel.background = element_rect(color = "black", fill = "transparent")) +
    labs(
      y = expression(paste("Absolute biomass (mg C ", m^{-3}, ")")),
      x = "Year",
      color = "Protist Maximum Body Size",
      shape = "Temperature Optima",
      fill = "Heatwave Timing:")



# Assemble vertical
f3_all <- (f3_a / f3_b / f3_c) + plot_layout(guides = "collect") & theme(legend.position = "right")
f3_all

Code

#  Futzing with the legend and arrangement

# # Pull out the legend separately if we want to reorganize
cop_leg <- as_ggplot(get_legend(f3_a))
prot_leg <- as_ggplot(get_legend(f3_c))
both_leg <- cop_leg / prot_leg  +  plot_layout(heights = c(3,1))

# # Rebuild
f3_rebuilt <- (f3_a | f3_c) / (f3_b | both_leg) & theme(legend.position = "none")
f3_rebuilt

Figure Change: No body size

Body Size isn’t talked about in the text so we could drop that and simplify the figures

Code

# Color and data for the heatwave boxes
hw_fill <- c("transparent", "#FDDBC770")


# Colors for thermal optima
manual_cols <- c(
  "16C" = gmri_cols("blue"),
  "20C" = gmri_cols("green"),
  "24C" = gmri_cols("yellow"),
  "28C" = gmri_cols("orange")
)



# Dominant Copepod
rect_df <- data.frame(
  hw_status = c("Year 0:\nPre-heatwave Equilibrium", "Year 1:\nSeasonal heatwave"),
  xmin = c(-0.5, 0.5),
  xmax = c(0.5, 1.5),
  ymin = rep(0, 2),
  ymax = rep(6, 2))
f3_a <-  ggplot(dom_new_ac) +
    geom_rect(
      data = rect_df,
      aes(xmin = xmin, xmax = xmax,
          ymin = ymin, ymax = ymax,
          fill = hw_status),
      color = "transparent",
      alpha = 0.9) +
    geom_line(
      aes(year, absolute_bio, group = taxa_id, color = temp_opt), 
      show.legend = F, 
      linewidth = 0.7, 
      alpha = 0.4) +
    geom_point(
      aes(year, absolute_bio, color = temp_opt), 
    size = 2) +
    scale_color_manual(values = manual_cols) +
    scale_fill_manual(values = hw_fill) +
    scale_x_continuous(
      limits = c(0, 7),
      expand = expansion(add = c(0.5,0.5)),
      breaks = c(0:9)) +
    guides(
      color = "none",
      fill = "none") +
    facet_grid(season~hw_season, scales = "free_y") +
    theme(
      legend.position = "right",
      panel.grid = element_blank(),
      panel.background = element_rect(color = "black", fill = "transparent"),
      strip.text.y = element_text(angle = 0)) +
    labs(
      title = "A. Active Copepods",
      y = expression(paste("Absolute biomass (mg C ", m^{-3}, ")")),
      x = NULL,
      color = "Thermal Optima",
      fill = "Heatwave Timing:")




# Passive Copepod
rect_df <- data.frame(
  hw_status = c("Year 0:\nPre-heatwave Equilibrium", "Year 1:\nSeasonal heatwave"),
  xmin = c(-0.5, 0.5),
  xmax = c(0.5, 1.5),
  ymin = rep(0, 2),
  ymax = rep(7, 2))
f3_b <- ggplot(dom_new_pc) +
    geom_rect(
      data = rect_df,
      aes(xmin = xmin, xmax = xmax,
          ymin = ymin, ymax = ymax,
          fill = hw_status),
      color = "transparent",
      alpha = 0.9) +
    geom_line(
      aes(year, absolute_bio, group = taxa_id, color = temp_opt), 
      show.legend = F, 
      linewidth = 0.7, 
      alpha = 0.4) +
    geom_point(
      aes(year, absolute_bio, color = temp_opt), 
    size = 2) +
    scale_color_manual(values = manual_cols) +
    scale_fill_manual(values = hw_fill) +
    scale_y_continuous(limits = c(0, 8) ) +
    scale_x_continuous(
      limits = c(0, 7),
      expand = expansion(add = c(0.5,0.5)),
      breaks = c(0:9)) +
    guides(
      shape = "none",
      color = "none",
      fill = "none") +
    facet_grid(season~hw_season, scales = "free_y") +
    theme(
      legend.position = "right",
      panel.grid = element_blank(),
      panel.background = element_rect(color = "black", fill = "transparent"),
      strip.text.y = element_text(angle = 0)) +
    labs(
      title = "B. Passive Copepods",
      y = expression(paste("Absolute biomass (mg C ", m^{-3}, ")")),
      x = NULL,
      color = "Thermal Optima",
      fill = "Heatwave Timing:")




# Protists
rect_df <- data.frame(
  hw_status = c("Year 0:\nPre-heatwave Equilibrium", "Year 1:\nSeasonal heatwave"),
  xmin = c(-0.5, 0.5),
  xmax = c(0.5, 1.5),
  ymin = rep(0, 2),
  ymax = rep(28, 2))
f3_c <-  ggplot(dom_new_pr) +
    geom_rect(
      data = rect_df,
      aes(xmin = xmin, xmax = xmax,
          ymin = ymin, ymax = ymax,
          fill = hw_status),
      color = "transparent",
      alpha = 0.9) +
    geom_line(
      aes(year, absolute_bio, group = taxa_id, color = temp_opt), 
      show.legend = F, 
      linewidth = 0.7, 
      alpha = 0.4) +
    geom_point(
      aes(year, absolute_bio, color = temp_opt), 
      size = 2) +
    scale_color_manual(values = manual_cols) +
    scale_fill_manual(values = hw_fill) +
    scale_x_continuous(
      limits = c(0, 7),
      expand = expansion(add = c(0.5,0.5)),
      breaks = c(0:9)) +
    guides(
      color = guide_legend(
        order = 1,
        title.position = "top", 
        title.hjust = 0, 
        direction = "vertical",
        ticks.colour = "black", 
        frame.colour = "black", 
        override.aes = list(size = rep(3, 4))),
      fill = guide_legend(
        order = 2,
        title.position = "top", 
        title.hjust = 0, nrow = 2,
        direction = "vertical",
        keyheight = unit(0.8, "cm"),
        override.aes = list(color = "black"))) +
    facet_grid(season~hw_season) +
    theme(
      legend.position = "right", 
      panel.grid = element_blank(),
      panel.background = element_rect(color = "black", fill = "transparent"),
      strip.text.y = element_text(angle = 0)) +
    labs(
      title = "C. Protists",
      y = expression(paste("Absolute biomass (mg C ", m^{-3}, ")")),
      x = "Year",
      color = "Thermal Optima",
      fill = "Heatwave Timing:")



# Assemble
f3_new <- (f3_a / f3_b / f3_c) + plot_layout(guides = "collect") & theme(legend.position = "right")
f3_new

Code

#  Futzing with the legend and arrangement

# # Pull out the legend separately if we want to reorganize
legend_gg <- as_ggplot(get_legend(f3_c) )

# # Rebuild
f3_nosize <- (f3_a | (f3_c+theme(legend.position = "none"))) / (f3_b | legend_gg) 
f3_nosize

Plots that follow in-text statements

Code

# about: the plots with 500 lines are so hard to
# make sense of. If the comparison is the different thermal preferences, let's highlight that


# What if we put everything into one table
# don't facet everything out
# use color for thermal preference

dom_all <- bind_rows(
  list(
    dom_new_ac,
    dom_new_pc,
    dom_new_pr)) %>% 
  mutate(taxa_id = str_c(functional_group, "-", taxa_id))

# Get the metadata
dom_all_meta <- distinct(dom_all, functional_group, taxa_id, temp_opt, min_size, max_size)

Statement 1: Community Composition - Thermal Opt

Overall, the plankton community is dominated by groups with temperature optima of 20 ˚C and 24 ˚C (Fig. 3).

Claires comment: this statement needs functional groups and should be annual summary.

Code

# C: Composition 

# Do an all season total
dom_all_totals <- bind_rows(
  list(
    # total up all the seasons within each hw season
    dom_all %>% 
      group_by(temp_opt, functional_group, hw_season) %>% 
      summarise(absolute_bio = sum(absolute_bio),
               .groups = "drop") %>% 
      mutate(season = "Overall Total"),
    # total up the seasons
    dom_all %>% 
      group_by(functional_group, temp_opt, hw_season, season) %>% 
      summarise(absolute_bio = sum(absolute_bio),
               .groups = "drop")
  )) %>% 
  mutate(season = factor(
    season, levels = c("Overall Total", "Spring", "Summer", "Autumn", "Winter"))) 




# table for hw event boxes
mhw_df <- data.frame(
  "x" = rep(1, 4), 
  "y" = rep(1, 4),
  season = c("Winter", "Autumn", "Spring", "Summer"),
  "hw_season" = str_c(c("Winter", "Autumn", "Spring", "Summer"), " Heatwave")) %>% 
  mutate(season = factor(
    season,
    levels = c("Overall Total", "Spring", "Summer", "Autumn", "Winter")))





# Make the plot
dom_all_totals %>% 
ggplot() +
    scale_fill_gmri() +
    geom_col(
      aes(
        y = functional_group, 
        x = absolute_bio, fill = temp_opt), 
      color = "gray40",
      alpha = 0.8, position = "fill") +
    scale_x_continuous(
      breaks = c(0.25, .5, 0.75), 
      labels = label_percent(),
      expand = expansion(add = c(0.05,0.05))) +
    facet_nested(hw_season ~ season ) +
    theme(strip.text.y = element_text(angle = 0),
          panel.grid = element_blank()) +
    labs(title = "Dominant Taxa Seasonal Composition",
         fill = "Temperature\nOptima",
         x = "Percent Abundance",
         y = "Functional Group")

Try doughnut plots, overall totals

Code

# Make the plot - doughnuts

dom_all_totals %>% 
  #filter(season == "Overall Total") %>% 
  ggplot() +
    scale_fill_gmri() +
  geom_arc_bar(
    aes(x0 = 0, y0 = 0, r0 = 0.2, r = 0.35, 
        fill = temp_opt, amount = absolute_bio), 
    stat = "pie") +
    scale_x_continuous(
      breaks = c(.5), 
      labels = label_percent(),
      expand = expansion(add = c(0.05,0.05))) +
    facet_nested(functional_group + season ~ hw_season ) +
    theme(
      strip.text.y = element_text(angle = 0),
      panel.grid = element_blank(),
      axis.line = element_blank(),
      axis.text = element_blank(),
      legend.position = "bottom", 
      legend.direction = "horizontal") +
    labs(title = "Dominant Taxa Seasonal Composition",
         fill = "Temperature\nOptima",
         x = "Biomass Percentage")

Statement 2: Community Composition - Thermal Opt

The temperature norms of dominant groups track the annual mean temperature, not the Sea Surface Temperature seasonality (Fig. 3),

This statement is not in original figure three at all, and speaks to whether or not species follow what the ambient temperature is or not.

Code

# C: Composition 

# Do an all season total
dom_all_yr_totals <- bind_rows(
  list(
    # total up all the seasons within each hw season
    dom_all %>% 
      group_by(year, temp_opt, hw_season) %>% 
      summarise(absolute_bio = sum(absolute_bio),
               .groups = "drop") %>% 
      mutate(season = "Overall Total"),
    # total up the seasons
    dom_all %>% 
      group_by(year, temp_opt, hw_season, season) %>% 
      summarise(absolute_bio = sum(absolute_bio),
               .groups = "drop")
  )) %>% 
  mutate(season = factor(
    season, levels = c("Overall Total", "Spring", "Summer", "Autumn", "Winter"))) 




# table for hw event boxes
mhw_df <- data.frame(
  "x" = rep(1, 4), 
  "y" = rep(1, 4),
  season = c("Winter", "Autumn", "Spring", "Summer"),
  "hw_season" = str_c(c("Winter", "Autumn", "Spring", "Summer"), " Heatwave")) %>% 
  mutate(season = factor(
    season,
    levels = c("Overall Total", "Spring", "Summer", "Autumn", "Winter")))





# Make the plot
dom_all_yr_totals %>% 
ggplot() +
    geom_col(
      aes(x = year, y = absolute_bio, fill = temp_opt), 
      alpha = 0.8, position = "fill", color = "gray40") +
    geom_col(
      data = mhw_df,
      aes(x,y), 
      color = "black", fill = "transparent", linewidth = 1) +
    scale_fill_gmri() +
    scale_y_continuous(breaks = c(0.25,.5, 0.75), 
                       labels = label_percent()) +
    scale_x_continuous(
      limits = c(-0.5, 6.5),
      expand = expansion(add = c(0.5,0.5)),
      breaks = c(0:6)) +
    facet_nested(hw_season ~ season ) +
    theme(strip.text.y = element_text(angle = 0)) +
    labs(title = "Dominant Taxa Seasonal Composition",
         fill = "Temperature\nOptima",
         x = "Year",
         y = "Abundance %")

Statement 3: Community Composition - Functional Groups

Examining community composition, heatwaves alter the order of dominant groups based on their relative contribution to the total biomass at the time (Fig. 3).

^This speaks to functional group change following MHW

If we combine the functional groups to one table we can do the community as a whole, and highlight the temperature optima of whatever and not focus so heavily on functional groups

Code

# So we're interested in composition


# Do an all season total
dom_fgroup_totals <- bind_rows(
  list(
    # total up all the seasons within each hw season
    dom_all %>% 
      group_by(year, functional_group, hw_season) %>% 
      summarise(absolute_bio = sum(absolute_bio),
               .groups = "drop") %>% 
      mutate(season = "Overall Total"),
    # total up the seasons
    dom_all %>% 
      group_by(year, functional_group, hw_season, season) %>% 
      summarise(absolute_bio = sum(absolute_bio),
               .groups = "drop")
  )) %>% 
  mutate(season = factor(
    season, levels = c("Overall Total", "Spring", "Summer", "Autumn", "Winter"))) 



# Make the plot
dom_fgroup_totals %>% 
  ggplot() +
  geom_col(
    aes(x = year, y = absolute_bio, fill =  functional_group), 
    alpha = 0.8, color = "gray40", position = "fill") +
  geom_col(
    data = mhw_df,
    aes(x,y), 
    color = "black", fill = "transparent", linewidth = 1) +
  scale_fill_gmri() +
  scale_y_continuous(
    breaks = c(0.25,.5, 0.75), 
    labels = label_percent()) +
  scale_x_continuous(
    limits = c(-0.5, 6.5),
    expand = expansion(add = c(0.5,0.5)),
    breaks = c(0:6)) +
  facet_nested(hw_season ~ season ) +
  theme(strip.text.y = element_text(angle = 0)) +
  labs(title = "Dominant Taxa Functional Group Composition",
       fill = "Functional Group",
       x = "Year",
       y = "Abundance %")

Change in year x from year x-1

Code

# B. Net-Change
dom_all %>% 
  group_by( year, temp_opt, hw_season, season) %>% 
  summarise(absolute_bio = sum(absolute_bio),
            .groups = "drop") %>%
  group_by(temp_opt, hw_season, season) %>% 
  arrange(year) %>% 
  mutate(bio_change = lag(absolute_bio, 1) - absolute_bio) %>% 
  ggplot() +
    geom_rect(
    data = rect_df[2,],
      aes(xmin = xmin, xmax = xmax,
          ymin = ymin, ymax = ymax),
      fill = hw_fill[2],
    color = "transparent", alpha = 0.9) +
    scale_fill_gmri() +
    geom_col(
      aes(x = year, y = bio_change, fill = temp_opt), 
      alpha = 0.8,
      color = "gray40") +
    facet_nested(hw_season ~ season) +
    theme(strip.text.y = element_text(angle = 0)) +
    labs(title = "B. Biomass Change")